Monte Carlo study of the magnetic critical properties 
of the two-dimensional Ising fluid 
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A two-dimensional fluid of hard spheres each having a spin ±1 and interacting via short-range 
Ising-like interaction is studied near the second order phase transition from the paramagnetic gas 
to the ferromagnetic gas phase. Monte Carlo simulation technique and the multiple histogram data 
analysis were used. By measuring the finite-size behaviour of several different thermodynamic quan- 
tities, we were able to locate the transition and estimate values of various static critical exponents. 
The values of exponents @/u and •y/u are close to the ones for the two-dimensional lattice Ising 
model. However, our result for the exponent v — 1.35 is very different from the one for the Ising 
universality class. 



CO 



I 

C 

O 

o 



> 

00 
On 



On 



-a 
c 

o 
o 



X 



I. INTRODUCTION 



Models with coupled translational and spin degrees 
of freedom attracted recently considerable attention be- 
cause they can describe several phenomena in amorphous 
ferromagnets 0, dilute magnetic alloys and dipolar liq- 
uids (III. Spins in such systems are not localized on 
lattice sites but are able to move. The systems with 
coupled spin and translational degrees of freedom ex- 
hibit a rich variety of phase transitions. Their phase 
diagrams were determined using the mean spherical ap- 
proximation, the mean-field theory, density functional 
methods and Monte Carlo (MC) simulation techniques 
However the critical behaviour and universality 
near phase transitions in these systems attracted only 
little attention. In ref. || the critical properties of the 
Heisenberg fluid near magnetic order-disorder transition 
were studied by MC simulations. The obtained critical 
exponents differ by a small but significant amount from 
the ones for the lattice Heisenberg model. The spin fluid 
systems resemble lattice-based spin models with annealed 
site dilution. The Blume-Capel model is an example of 
the lattice-based Ising model with an annealed site dilu- 
tion p]. The density of annealed sites in this model is, 
however, not fixed but can fluctuate around an average. 

The phase diagram of the 2d Ising fluid studied in 
this paper obtained within the mean-field approxima- 
tion described in the ref. is shown in figure [|. For 
high temperatures and low densities there exists critical 
line separating a paramagnetic gas phase from a ferro- 
magnetic gas phase. The critical line finishes at lower 
temperatures at the tricritical point. The similar phase 
diagram was obtained in the 2d Blume-Capel model. 
Recently MC simulations were used to investigate the 



tricritical point properties of a 2d Ising fluid and 2d 
Blume-Capel model 0. It was shown that both models 
belong to the same tricritical universality class. 
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FIG. 1. Mean-field phase diagram for the Ising fluid with 
the coupling constant K and the density p. The line separat- 
ing the paramagnetic and ferromagnetic gas phases is given 
by K — 1/ (2np) . The point indicates the location of the 
phase transition for p = 0.4 obtained in this paper from MC 
simulations. 
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The aim of this paper is to present the results obtained 
from MC simulations of the 2d spin fluid with short-range 
Ising-like interactions near the second-order phase transi- 
tion from the paramagnetic gas to the ferromagnetic gas 
phase far from the tricritical point. We performed simu- 
lations in systems of different sizes at a constant particle 
density and for four different temperatures. We analyzed 
our data combining the multiple histogram technique ||] 
with finite-size scaling (FSS) ||[l0| to obtain estimates 
for the critical temperature and exponents. In the Sec. 
II we describe the Ising fluid model and technical aspects 
of the simulations. Determination of critical exponents 
and the location of phase transition are given in Sees. Ill 
and IV, while Scc.V summarizes our results. 
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FIG. 2. The magnetization versus the coupling constant 
K for the five system sizes. The number of particles in the 
system is indicated. The points are averages over individual 
simulations. The curves results from the multiple histogram 
technique. Error bars are smaller than the symbol size. 



II. THE MODEL AND SIMULATION DETAILS 

We consider a system which consists of particles of di- 
ameter a in two-dimensional space. The internal degrees 
of freedom of each particle are described by an Ising spin 
and there is an exchange coupling between spins given 
by the Yukawa interaction. The system Hamiltonian is: 
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where the interaction potential has the following form, 

M \ - / 00 if r < cr , , 

W - 1 -K 2 exp(-{r-a)/a) if r > a [Z> 



The parameter K is the ratio of the coupling energy 
to the thermal energy. In this paper we consider the 
ferromagnetic case (K > 0). 




FIG. 3. The magnetic susceptibility versus the coupling 
constant K for the five system sizes. The number of par- 
ticles in the system is indicated. The points are averages over 
individual simulations. The curves results from the multiple 
histogram technique. Error bars are omitted when smaller 
than the symbol size. 

The MC simulations were performed at a constant par- 
ticle density p — 0.4. We studied five systems with the 
number of particles N equal to 128, 256, 512, 768 and 
1024. The periodic boundary conditions and the mini- 
mum image convention were applied during simulations 
tyj. The interaction potential <p(r) was cut at a dis- 
tance 6.3246tr. This value was chosen in order to divide 
the simulation cell into 16 sub-cells for the system with 
256 particles. In order to speed up simulations we used 
the method of linked lists of neighbours tyj. We ap- 
plied the same simulation algorithm as it is described 
in ref. In the present work we have not included 
any long-range correction to the cuttof procedure as it 
was done in |J. The maximum position displacement of 
particles was chosen in such a way, that the acceptance 
ratio of the trial moves was around 0.5. The number 
of MCS/N (Monte Carlo steps per particle) discarded 
at the beginning of the simulation was larger than 10 4 . 
For each system size the simulations were performed for 
four values of the parameter K. These values were the 
following: K = 0.48, 0.5, 0.52 and 0.54 for a system with 
128 particles, K = 0.48, 0.515, 0.53 and 0.54 for a system 
with 256 particles and K = 0.485,0.515,0.53 and 0.54 
for systems with 512,768 and 1024 particles. These val- 
ues of K were chosen, because the range of temperatures 
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where reliable extrapolation of the behaviour of physical 
quantities could be performed included positions of the 
maxima of specific heat and magnetic susceptibility and 
the position of the phase transition in the bulk system. 
The data were stored at intervals of 10 MCS/N and the 
total number of updates was 2 • 10 6 MCS/N for every K 
value. 

Let's denote by M the magnetization per parti- 
cle of the system defined as M = (%2iSi)/N. Wc 
study critical behaviour of the following quantities 
jn^ , [l3| : the mean absolute value of the magnetiza- 
tion < \M\ >, the magnetic susceptibility defined as 
X = KN{< M 2 > - < \M\ > 2 ) the fourth-order magne- 
tization cumulant defined as U = 1— < Af 4 > /(3 < 
M 2 >), the mean energy of the system < H > and 
the specific heat C v — (< H 2 > - < H > 2 )/N, where 
< ... > denotes canonical ensemble average. We also con- 
sider the quantities like the derivatives d In < \M\ > JdK 
and 9 In < M 2 > jdK. In the Heisenberg fluid || and 
3d lattice Ising model |ll| these derivatives were used to 
extract the critical exponent v. 
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FIG. 4. The fourth-order magnetization cumulant versus 
the coupling constant K for the five system sizes. The num- 
ber of particles in the system is indicated. The points are 
averages over individual simulations. The curves results from 
the multiple histogram technique. Error bars are omitted 
when smaller than the symbol size. 

The thermodynamic properties of a system, in a wide 
temperature range, can be obtained by performing sev- 
eral MC simulations at different temperatures and com- 
bining them by the application of the multiple histogram 
technique ||. This technique allows reliable extrapola- 
tion of MC results to the values of K where the interesting 
positions of the maxima shown by some of the quantities 



defined above are located. In figures g, |3j and |4| we show 
the dependence of quantities < \M\ >, \ and U on K. 
The results obtained directly from MC simulations per- 
formed at selected K values are also shown in these figures 
by points. The error bars were obtained by calculating 
block averages of 10 5 MCS/N data points and computing 
the standard deviation from these block averages. One 
can notice that the multiple histogram extrapolations are 
within calculated error bars. 
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FIG. 5. The plot of the maxima of *M<jm> an d 81n <* f2 > 
versus the linear system size L. Error bars are smaller than 
the symbol size. The straight lines are fits to the data. Their 
slope is 1/v. The obtained values of the exponent v are indi- 
cated. 



III. ESTIMATES OF THE EXPONENT v AND 
THE LOCATION OF THE PHASE TRANSITION 
IN THE BULK SYSTEM 

The critical exponent v characterizing the divergence 
of the correlation length near the second-order phase 
transition [|l0| can be extracted by considering the scal- 
ing behaviour of the derivatives 9 In < \M\ > jdK and 
d In < M 2 > j dK [^|Jl^1 . These derivatives can easily be 
computed using the following identity: 

8 < I M n I > 1 

' = (< \M n \H >-< \M n \ >< H >) 
oK K 

(3) 

Let's denote by L the lengt h of o ne side of the simu- 
lation box. In our case L = \/N/ p. The dependence of 
these derivatives on K has the maximum which should 
scale with the system size as L 1 ^ . This method of es- 
timation of the exponent v is very convenient, because 
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it can be done without any consideration of the criti- 
cal coupling K c in the bulk system. We show in the 
figure ^| maximum values of d In < \M\ > /dK and 
9 In < M 2 > /dK in systems of different sizes together 
with the fitted straight lines. The goodness of fits Q 
jf4| are Q = 0.36 and Q = 0.32 for the derivative of 
In < | M | > and In < M 2 > respectively. The slope 
of these lines provides estimates for v, and we obtained 
v = 1.36 ±0.02 for the maxima of d In < \M\ > /dK and 
v = 1.34 ± 0.02 for the maxima of d In < M 2 > /dK. 
These values are higher than the value v = 1 in the two- 
dimensional lattice Ising model. 
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FIG. 6. Size dependence of the location of maxima K C (L) 
for several quantities indicated in the figure in finite systems 
containing 128, 256, 512, 768 and 1024 particles. The straight 
lines are fits to Eq. ^ with the exponent v = 1.35. The esti- 
mated values of the critical coupling K c in the bulk system 
obtained from these fits are shown for each quantity. 

The critical coupling K c in the bulk system is usually 
determined using the Binder fourth-order magnetization 
cumulant crossing technique |5|,|l^,[l3| . Finite size scaling 
predicts that for sufficiently large systems if we plot U 
versus K for different choices of L, these curves should 
have a unique intersection point U* . The value of K 
where this occurs is the value of the critical coupling K c . 
This value is not biased by any assumptions about critical 
exponents. For smaller systems there are corrections to 
FSS and the intersection point between any two curves 
U vs. K corresponding to systems with side lengths L 
and L' depends on L and L' Jig ]. In Table I we give 
the coordinates of intersection points for different pairs 
of systems. The values in the table have statistical er- 
rors larger than the expected correction terms to FSS. 
Because of this and a small number of systems studied, 
we were not able to extract K c by the extrapolation pro- 



cedures given in ref. [[l3|.The critical coupling K c and the 
intersection value U* we calculated as the average from 
the values in the table. We obtained K c = 0.535 ± 0.002 
and U* = 0.618 ± 0.003. We have excluded from the av- 
erage the last three rows of the table because for systems 
with small difference in sizes even a small shift in the cu- 
mulant lines can produce a considerable error in the co- 
ordinates of intersection points. The estimated common 
value of the cumulant is only slightly larger than the value 
U* = 0.611 ±0.001 obtained for the two-dimensional lat- 
tice Ising model Q. 

The value of K c can also be determined from the size- 
dependent shifting of the peak of different thermody- 
namic quantities. In finite systems the quantities like 
e.g. the specific heat CV, the magnetic susceptibil- 
ity x, d < \M\ > /dK , dhx < \M\ > /dK and 
9 In < M 2 > /dK, have maxima as a function of K 
[p|, p~0|Jl3| . The location of the maximum K C (L,A) de- 
pends both on the system size L and on the quantity 
A considered. FSS predicts the following dependence of 
K C (L, A) on the system size |p|,^3|: 

K C {L,A) =K c + aL- 1 / u (4) 

with the omitted corrections to FSS. The constant a de- 
pends in magnitude and sign on the particular quantity 
considered. In order to determine K c from this equation 
it is necessary to have both an accurate estimate of the 
exponent v and accurate values of K C {L,A). In figure 
^ we plot estimates of K C (L,A) for different quantities 
as a function of L~ x l v . The lines in this figure were ob- 
tained by the least square fits of the data to Eq. || with 
v = 1.35, the average of previously determined values. 
One can notice that values of K c obtained from the fit 
agree well with the estimated value K c = 0.535 estimated 
above from the cumulant intersection points. The value 
of the exponent v different than 1.35 will lead, of course, 
to different estimates of K c which are inconsistent with 
the value extracted from cumulant intersection points. 



IV. ESTIMATES OF EXPONENT RATIOS (3/v 

AND y/v 



The exponent ratio (5/v can be obtained from the 
finite-size scaling behaviour of \M either at K = K c or at 
the value of K where the derivative d < \M\ > /dK has 
the maximum. FSS predicts that \M\ at these K values 
should obey the relation M ~ L~ l3 l v . Figure shows 
the plots corresponding to this relation. The straight 
lines in this figure were obtained using the least-square 
fitting routine. The exponent j3 /v was determined from 
the slope of the lines. We obtained: /3/u = 0.141 ± 0.005 
(Q = 0.74) at K c and (3/u = 0.120 ± 0.006 (Q = 0.77) at 
the maximum of d < \M\ > /dK. 
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FIG. 7. The plot of the magnetization at the estimated 
value of the critical coupling in the bulk system K c = 0.535 
and at the maximum of d< g^ > , versus the linear system size 
L. The straight lines are fits to the data. Their slope is 
—/3/v. The obtained values of the ratio of exponents fi/v are 
indicated. 
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FIG. 8. The plot of the magnetic susceptibility at the es- 
timated value of the critical coupling in the bulk system 
K c — 0.535 and at the maximum, versus the linear system 
size L. The straight lines are fits to the data. Their slope is 
7/V. The obtained values of the ratio of exponents 7/v are 
indicated. 



The exponent ratio j/v can be determined from the 



finite-size scaling behaviour of the maximum of the mag- 
netic susceptibility x max and of the value x(-^c) of the 
magnetic susceptibility at K — K c . According to FSS 
these quantities are expected to vary with system size like 
ipl v , Figure || displays the finite-size scaling behaviour 
of Xmax and x(-Kc)- We estimated the value of the expo- 
nent ratio 7/zv from slopes of fitted straight lines. From 
values Xmax we obtained 7/1/ = 1.74 ± 0.02 (Q = 0.31), 
whereas from values x(Kc) we obtained 7/1/ = 1.76±0.05 
(Q = 0.88). 



V. CONCLUSIONS 

We have studied the critical behaviour of the two- 
dimensional Ising fluid at a density p — 0.4 near the 
second-order phase transition from the paramagnetic gas 
to the ferromagnetic gas phase. The multiple histogram 
technique was applied in order to combine data obtained 
from four different MC simulations. The largest sys- 
tem we studied consisted of 1024 particles. The critical 
exponent v was determined from FSS behaviour of the 
derivatives 9 In < |M| > jdK and 9 In < M 2 > jdK. 
We obtained v ps 1.35. This value is much higher than 
the value v = 1 obtained for the two-dimensional Ising 
lattice model. The critical coupling K c in the bulk sys- 
tem was obtained by considering size-dependent shifting 
of the maxima of several quantities using the estimated 
value of the exponent v. This value was found to be 
consistent with the location of the fourth-order magne- 
tization cumulant crossing points i.e. with K c « 0.535. 
The estimated common value of the cumulant at K = K c 
is U* = 0.618 what is slightly higher than the value 
U* = 0.611 obtained in the two-dimensional lattice Ising 
model. All the cumulant intersection values given in table 
I are above the value 0.611 except the value determined 
for the pair of systems [768, 1024]. As the statistical er- 
ror increases with the system size and considerable error 
is expected in the determination of cumulant intersec- 
tion points for pairs of systems with small difference in 
sizes we have neglected the last three rows of table I in 
estimating K c and U* . The ratio of critical exponents 
j3 jv and j/v we determined from the FSS behaviour of 
the magnetization and magnetic susceptibility. They are 
consistent with values obtained in the two-dimensional 
lattice Ising model. In annealed diluted magnets the ex- 
ponents are renormalized if the pure specific-heat expo- 
nent a is positive and do not change if it is negative 
jnj. As for the pure lattice Ising model a — 0, we ex- 
pect for the annealed diluted Ising models the pure Ising 
exponents. Our results disagree with this prediction and 
suggests a renormalization of exponents which leaves (3/v 
and 7/1/ unchanged and changes v. We were unable to 
consider large system sizes with good statistics to deter- 
mine possible corrections to finite-size scaling. However, 
the presented self-consistency in the determination of K c 
using the estimated exponent v convinces one that the 
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obtained values of v and K c are correct. 

A behaviour similar to ours was reported in studies 
of the two-dimensional quenched site diluted Ising model 
fyl. These studies suggest that the v exponent and the 
Binder cumulant value, U* increases with the degree of 
disorder. These was interpreted as the verification of the 
weak universality scenario (l9) also seen to apply to other 
models [fL8[ . We plan to study our model with different 
densities to see if the measured exponents approach the 
pure Ising values with increasing density. 

TABLE I. The values of fourth-order magnetization cumu- 
lant U C ross(Ni, N2) and the coupling K cr oss(Ni, N2) at the 
intersection point for different pairs of systems [Ni,N2] hav- 
ing Ni and N2 particles. 



[Ni,N 2 ] 


K cross {N ly N 2 ) 


U C ross{Nl, 


[128, 256] 


0.5327 


0.6153 


[128, 512] 


0.5331 


0.6158 


[128, 768] 


0.5364 


0.6193 


[128, 1024] 


0.5358 


0.6186 


[256, 512] 


0.5334 


0.6164 


[256, 768] 


0.5380 


0.6224 


[256, 1024] 


0.5368 


0.6209 


[512, 768] 


0.5449 


0.6336 


[512, 1024] 


0.5398 


0.6269 


[768, 1024] 


0.5324 


0.6105 



N 2 ) 
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